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ABSTRACT 

A finite element approach using NASTRAN is developed for solving time-dependent fluid- 
structure interaction problems, with emphasis on the transient scattering of acoustic \va\es from 
submerged elastic structures. Finite elements are used for modeling both structure and fluid domains 
to facilitate the graphical display of the wave motion through both media. For the fluid, the use of 
velocity potential as the fundamental unknown results in a symmetric matrix equation. The approach 
is illustrated for the problem of transient scattering from a submerged elastic spherical shell subjected 
to an incident tone burst. The use of an analogy between the equations of elasticity and the wave 
equation of acoustics, a necessary ingredient to the procedure, is summarized. 


INTRODUCTION 

Computational structural acoustics is concerned with the prediction of the acoustic pressure 
field radiated or scattered by submerged structures subjected to either mechanical or external (fluid) 
excitation. When the excitation is time-harmonic, the most common numerical approach for solving 
the interaction problem is to couple a finite element model of the structure with a boundary element 
model of the surrounding fluid (Ref. 1-8). Other fluid modeling approaches have included finite 
element (Ref. 9-20), combined finite element/analytical (Ref. 21-23), and T-matrix (Ref. 24-26). 

For time domain (transient) analysis, there are several computational approaches which can be 

used: 

• the transformation of frequency domain results to the time domain using the Fourier transform 

• the use of a fluid loading approximation such as the doubly asymptotic approximation (DAA) 
(Ref. 27) 

• the time domain boundary element approach, w'hich models the fluid with the retarded potential 
integral equation (Ref. 28-31) 

• the fluid finite element approach, which models the exterior fluid domain with finite elements 
truncated at a finite distance from the structure and terminated with an approximate radiation 
boundary condition to absorb outgoing weaves (Ref. 9-20) 

To our knowledge, the retarded potential integral equation has been used only for special geometries 
(e.g., axisymmetry) because of the method’s relatively high computational cost. The DAA approach, 
which has been used successfully in underwater shock analysis (Ref. 32-34), may not be adequate for 
transient acoustics, where the interest is in the response in the fluid as well as in the structure. The 
principal computational trade-off between the fluid finite element approach and the other three 
approaches is that the finite element approach yields large, banded matrices, w^hereas the other three 
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approaches (which depend on boundary element calculations) yield smaller, densely-populated 
matrices. This trade-off often favors the finite element approach for long, slender structures like 
ships which are "naturally banded." In addition, of the four approaches listed, only the fluid finite 
clement approach has directly available an explicit fluid mesh which can be used for graphical display 
of the wave motion through the fluid. Since a significant part of our interest involves the display of 
wave propagation through both structure and fluid, we therefore formulate the transient acoustics 
problem using the fluid finite element approach. The principal drawbacks to a fluid finite element 
approach are the need for an approximate radiation boundary condition at the outer fluid boundary, 
the requirements on mesh size and extent (sometimes leading to frequency-dependent fluid meshes 
(Ref. 17)), and the difficulty of generating the fluid mesh. 

Dynamics problems involving the interaction between an elastic structure and an acoustic fluid 
have been formulated for finite element solution using pressure (Ref. 9,10), fluid particle 
displacement (Ref. 11-13,15,17), displacement potential (Ref. 16), and velocity potential (Ref. 18,19) 
as the fundamental unknown in the fluid region. In three dimensions, the pressure and displacement 
formulations result in, respectively, one and three degrees of freedom per finite element mesh point. 
Thus the pressure approach has the advantage of fewer unknowns and a smaller overall matrix profile 
or bandwidth if the grid points are properly sequenced. On the other hand, the displacement 
approach results in symmetric coefficient matrices (in contrast to the pressure formulation, for which 
the matrices are nonsymmetric) and a fluid-structure interface condition which is easier to implement 
with general purpose finite element computer programs. However, the displacement approach also 
suffers from the presence of spurious resonances (Ref. 15), a situation which can be bothersome in 
time-harmonic problems, either forced or unforced. The principal disadvantage of the pressure 
formulation, nonsymmetric coefficient matrices, can be removed merely by reformulating the pressure 
solution approach so that a velocity potential rather than pressure is used as the fundamental 
unknown in the fluid region (Ref. 18). For some situations, particularly time-harmonic problems 
involving damped systems and time-dependent problems, significant computational advantages result. 

The principal goal of this paper is to develop in detail the symmetric velocity potential 
formulation for application to the specific problem of transient acoustic scattering from submerged 
elastic structures. Previously (Ref. 18), the symmetric potential formulation was described only in 
general terms for a wide class of fluid-structure interaction problems with no details concerning 
specific types of applications such as vibrations, shock response, or acoustic scattering. 

From an engineering point of view, it is convenient to be able to make use of existing general 
purpose finite element codes such as NASTRAN, because of their wide availablity, versatility, 
reliability, consultative support, and abundance of pre- and postprocessors. Thus the next section 
summarizes an analogy between the equations of elasticity and the wave equation of acoustics. Such 
an analogy allows the coupled structural acoustic problem to be solved with standard finite element 
codes. 


STRUCTURAL-ACOUSTIC ANALOGY 

Since we wish to solve the coupled structural acoustic problem using standard finite element 
codes, we summarize here the application of such codes to the wave equation of acoustics (Ref. 
35,36), 

V 2 p = p/c 2 , W 

where V 2 is the Laplacian operator, p is the dynamic fluid pressure, c is the wave speed, and dots 
denote partial differentiation with respect to time. 

On the other hand, the x-component of the Navier equations of elasticity, which are the 
equations solved by structural analysis computer programs, is 
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( 2 ) 


_X+2G , X+G , \ 1 r P •• 

q u, xx + U >yy + u >zz + q ( v >xy + w >xz) + ^"fx = ^ u ’ 

where u, v, and w are the Cartesian components of displacement, X is a Lame' elastic constant, G is 
the shear modulus, f x is the x-component of body force per unit volume (e.g., gravity), p is the mass 
density, and commas denote partial differentiation. 

A comparison of Eqs. 1 and 2 indicates that elastic finite elements can be used to model scalar 
pressure fields if we let u, the x-component of displacement, represent p, set v = w = 0 everywhere, 

fx = 0, and X = -G. For three-dimensional analysis, the engineering constants consistent with this last 
requirement are (Ref. 36) 

E e = 10 20 G e , p e = G e /c 2 , ( 3 ) 

where the element shear modulus G e can be selected arbitrarily. The subscript "e" has been added to 
these constants to emphasize that they are merely numbers assigned to the elements. 

A variety of boundary conditions may also be imposed. At a pressure-release boundary, p = 0 
is enforced explicitly like other displacement boundary conditions. For gradient conditions, the 
pressure gradient dp/dn is enforced at a boundary point by applying a "force” to the unconstrained 
DOF at that point equal to G e Adp/5n, where A is the area assigned to the point and n is the outward 
normal from the fluid region (Ref. 36). For example, the plane wave absorbing boundary condition 



is enforced by applying to each point on the outer fluid boundary a "force” given by -(G e A/c)p. 

Since this force is proportional to the first time derivative of the fundamental solution variable p, 
this boundary condition is imposed in the analogy by attaching to the fluid DOF a "dashpot" of 
constant G e A/c. The Neumann condition dp/dn = 0 is the natural boundary condition under this 
analogy. The next higher order local radiation boundary condition, the curved wave absorbing 
boundary condition (Ref. 20,37) 

in = _i _ £. 

dn or’ ( 5 ) 


where r is the radius of the boundary, is enforced under the analogy by attaching in parallel both a 
"dashpot" and a "spring" between each boundary point and ground. 

At a fluid-structure interface (an accelerating boundary), momentum and continuity 
considerations require that 


in = 

dn 




( 6 ) 


where n is the normal at the interface, p is the mass density of the fluid, and ii n is the normal 
component of fluid particle acceleration. Under the analogy, this condition is enforced by applying to 
the fluid DOF at a fluid-structure interface a "force" given by — (G e />A)ii n . 

To summarize, the wave equation, Eq. 1, can be solved with elastic finite elements if the three- 
dimensional region is modeled with 3-D solid finite elements having material properties given by Eq. 

3, and only one of the three Cartesian components of displacement is retained to represent the scalar 
variable p. In Cartesian coordinates, any of the three components can be used. The solution of 
axisymmetric problems in cylindrical coordinates follows the same approach except that the z- 
component of displacement is the only one which can be used to represent p (Ref. 36). 
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SCATTERING FROM ELASTIC STRUCTURES 

In the scattering problem, a submerged elastic body is subjected to a plane wave incident 
loading, as shown in Fig. 1. For the time-harmonic case, the excitation has a single circular 
frequency w. For the time-dependent (transient) case of interest here, the pi escribed pressure 
loading is an arbitrary function of time. Without loss of generality, we can assume that the incident 
wave propagates in the negative z direction. The speed of such propagation is c, the speed of sound 
in the fluid. 



Fig. 1. The scattering problem. 

Within the fluid region, the total pressure p satisfies the wave equation, Eq. 1. Since the 
incident free-field pressure p; is known, it is convenient to decompose the total pressure p into the 
sum of incident and scattered pressures 

( 7 ) 

p = Pi + Ps> 

each of which satisfies the wave equation. (By definition, the incident free-field pressure is that 
pressure which would occur in the fluid in the absence of any scatterer.) 

We now formulate the problem for finite element solution. Consider an arbitrary, submerged, 
three-dimensional elastic structure subjected to either internal time-dependent loads or an external 
time-dependent incident pressure. If the structure is modeled with finite elements, the resulting 
matrix equation of motion for the structural degrees of freedom (DOF) is 

Mii + Bu + Ku = F — GAp, 

where M, B, and K are the structural mass, viscous damping, and stiffness matrices (dimension s x 
s), respectively, u is the displacement vector for all structural DOF (wet and dry) in terms of the 
coordinate systems selected by the user (s x r), F is the vector of applied mechanical forces applied 
to the structure (s x r), G is the rectangular transformation matrix of direction cosines to transform a 
vector of outward normal forces at the wet points to a vector of forces at all points in the coordinate 
systems selected by the user (s x f), A is the diagonal area matrix for the wet surface (f x f), p is the 
vector of total fluid pressures (incident + scattered) applied at the wet grid points (f x r), and dots 
denote differentiation with respect to time. The pressure p is assumed positive in compression. In 
the above dimensions, s denotes the total number of independent structural DOF (wet and dry), f 
denotes the number of fluid DOF (the number of wet points), and r denotes the number of load 
cases. If first order finite elements are used for the surface discretization, surface areas, normals, 
and the transformation matrix G can be obtained from the calculation of the load vector resulting 
from an outwardly directed static unit pressure load on the structure’s wet surface. The matrix 
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product GA can then be interpreted as the matrix which converts a vector of negative fluid pressures 
to structural loads in the global coordinate system. The last two equations can be combined to yield 

Mu' + Bu + Ku + G Ap s = F - GAp s . (9) 

A finite element model of the fluid region (with scattered pressure p s as the unknown) results in 
a matrix equation of the form 

QPs + Cp s + Hp s = F<p\ (10) 

where p s is the vector of scattered fluid pressures at the grid points of the fluid region, Q and H are 
the fluid "inertia" and "stiffness" matrices (analogous to M and K for structures), C is the "damping" 
matrix arising from the radiation boundary condition (Eq. 4), and F (p) is the "loading" applied to fluid 
DOF due to the fluid-structure interface condition, Eq. 6. Using the analogy described in the 
preceding section, structural finite elements can be used to model both structural and fluid regions. 
Material constants assigned to the elastic elements used to model the fluid are specified according to 
Eq. 3. In three dimensions, elastic solid elements are used (e.g., isoparametric bricks for general 3-D 
analysis or solids of revolution for axisymmetric analysis). 


At the fluid-structure interface, Eqs. 6 and 7 can be combined to yield 

3p s 

- f(n -' - “*>• 


(in 


where n is the outward unit normal, and ii ni and u„ are, respectively, the incident and total outward 
normal components of fluid particle acceleration at the interface. Thus, from the analogy, we impose 
the fluid-structure interface condition by applying a "load” to each interface fluid point given by 

F (p) = -pG e A (u ni - u' n ), (12) 

where the first minus sign is introduced since, in the coupled problem, we choose n as the outward 
normal from the structure into the fluid, making n an inward normal for the fluid region. The normal 
displacements u n are related to the total displacements u by the same rectangular transformation 
matrix G used above: 


u n = G T u, (13) 

where the superscript T denotes the matrix transpose. Eqs. 10, 12, and 13 can be combined to yield 

QPs + Cp s + Hp s - pG e (GA) T ii = -pG e Aii ni . (14) 

Since the fluid-structure coupling terms in Eqs. 9 and 14 are nonsymmetric, we symmetrize the 
problem (Ref. 18) by using a new fluid unknown q such that 


q = /p s dt, q = Ps- 
o 

If Eq. 14 is integrated in time, and the fluid element "shear modulus" G e is chosen as 
G e = -1/p, 

the overall matrix system describing the coupled problem can be written as 


M 

0 



B 

(GA) T 


(GA) 

C 



+ 



Jf — (GA) P jl 
|Av ni j> 


(15) 

(16) 
(17) 


where v ni (=u ni ) is the outward normal component of incident fluid particle velocity. 

The new variable q is, except for a multiplicative constant, the velocity potential <j>, since 
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(18) 


P = -p4>- 

Eq. 17 could also be recast in terms of <j> rather than q as the fundamental fluid unknown, but no 
particular advantage would result. In fact, the use of q rather than 4> has the practical advantage that 
the fluid pressure can be recovered directly from the finite element program as the time derivative 
(velocity) of the unknown q. 

To summarize, both structural and fluid regions are modeled with finite elements. For the fluid 
region, the material constants assigned to the finite elements are 

E e = —10 20 /p, G e = -1 Ip, = unspecified, p e = -l/(/?c 2 ), (19) 

where E e , G e , and Pc are the Young’s modulus, shear modulus, Poisson’s ratio, and mass density, 
respectively, assigned to the fluid finite elements. The properties p and c above are the actual density 
and sound speed for the fluid medium. The radiation boundary condition used is the plane wave 
approximation, Eq. 4, which appears to be adequate if the outer fluid boundary is sufficiently far from 
the structure (Ref. 17). With this boundary condition, matrix C in Eq. 17 arises from dashpots 
applied at the outer fluid boundary with damping constant -A/(pc) at each grid point to which the 
area A has been assigned. At the fluid-structure interface, matrix GA is entered using the areas (or 
areal direction cosines) assigned to each wet degree of freedom. (Recall that GA can be interpreted 
as the matrix which converts a vector of negative fluid pressures to structural loads in the global 
coordinate system.) 

The right-hand side of Eq. 17 can be simplified further since, for plane waves propagating in the 
negative z direction at speed c, the incident free-field pressure and incident fluid particle velocity in 
the z direction are related by (Ref. 38) 


Then, like in Fig. 1, if we define 8 as the angle between the normal n and the positive z axis, 

v ni = v zi cos0 = — picos 9 /{pc). (21) 

For plane waves, the z component of the free-field fluid particle velocity v zi is the same at all points in 
space except for a time delay, which depends only on the z coordinate of the points. 


Thus, Eq. 17 can alternatively be written 


M O' 

ful 


B (GA) 

w + 

K 0 

u _ 

.0 Q. 

W 

‘ + 

(GA) T C 

w 

0 H 

w 


F - (GA)pi 
Ap; cos 8 /{pc) I 


( 22 ) 


This is the form of the equations which we will use to solve the transient scattering problem. The 
right-hand side, which has nonzero contributions for both structure and fluid interface points, 
depends only on the incident free-field pressure at the fluid-structure interface. For scattering 
problems, the mechanical load F is zero. For radiation problems, F is nonzero, and the incident 
pressure p; vanishes. 


We note that the structural and fluid unknowns are not sequenced as perhaps implied by the 
partitioned form of Eq. 22. The coupling matrix GA is quite sparse and has nonzeros only for matrix 
rows associated with the structural DOF at the fluid-structure interface and columns associated with 
the coincident fluid points. Thus, the grid points should be sequenced for minimum matrix 
bandwidth or profile as if the structural and fluid meshes comprised a single large mesh. As a result, 
the structural and fluid grid points will, in general, be interspersed in their numbering, and the system 
matrices will be sparse and banded. 
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EXAMPLE: SCATTERING FROM A SUBMERGED SPHERICAL SHELL 

The validation of the procedure described above was made by comparing the finite element 
prediction of the time history of the structural response of a spherical shell subjected to a step 
incident pressure loading with the series solution (Ref. 28,39). These results will not be presented 
here. Instead, we will illustrate the approach by calculating the transient response of a submerged, 
thin-walled, evacuated spherical shell subjected to a brief tone burst, as illustrated in Fig. 2. For 
convenience, we nondimensionalize lengths to the shell mean radius a, velocities to the fluid sound 
speed c, and pressures to the fluid bulk modulus pc 2 . Thus, nondimensional time becomes ct/a. The 
particular problem solved was a 2% thick steel shell immersed in water. Hence, in nondimensional 
units, the shell properties are thickness = 0.02, Young’s modulus = 96.9, Poisson’s ratio = 0.3, and 
density = 7.79. 


VACUUM 



Fig. 2. Scattering from a submerged spherical shell. 


The incident free-field pressure Pj(z,t) is given by 


Pi(x,y,z,t) = p;(t + — ), 
c 


where (Fig. 3) 
Pi(t) = 


Po (1 — cos cot) /2, 

— po COS U>t, 

-Po (1 + cos wt)/2, 

0 , 


0 < wt < 7T 
7T < Wt < (n-l)7T 
(n— l)7r < wt < n7r 
otherwise. 


(n odd) 


(23) 


(24) 


For this problem, p Q =l, n=5, and wa/c= 7 r. 

Since this problem is axisymmetric, it was modeled for finite element solution using 
NASTRAN’s conical shell elements (CONEAX) for the shell and triangular ring elements 
(TRIAAX) for the fluid. A typical fluid mesh is shown in Fig. 4, where the shell is coincident with 
the inner semi-circle of fluid grid points. The actual mesh used to generate the results which follow 
had the outer fluid boundary located at 8 radii, used 24 elements along the inner radius between the 
poles and 56 elements in the radial direction, resulting in a total of 25 structural grid points, 6213 fluid 
grid points, 24 CONEAX elements, 12096 TRIAAX elements, and 6288 independent degrees of 
freedom. For the direct time integration, 800 nondimensional time steps of size 0.025 were used. 

Results are presented for both velocity response of the shell and scattered pressure response in 
the fluid. Fig. 5 shows plots of time histories of shell velocity in the z direction for the point first 
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Fig. 3. The incident pressure pj(t) (Eq. 24). 



Fig. 4. Typical finite element mesh. 

impacted by the pressure wave (0=0) and the back side pole (0=180 degrees). We observe from the 
figure a significant oscillation in the back side of the shell. Back-scattered pressure time histories are 
displayed in Fig. 6 at 3 and 5 radii from the origin. As expected, the scattered pressure at fluid 
points is zero until the wave has had time to travel (at unit nondimensional speed) from the spherical 
shell. Since the two points displayed are located 2 and 4 radii from the shell, the nondimensional 
time delays for the scattered pressure wave to arrive are 2 and 4, respectively. 

DISCUSSION 

A practical procedure has been presented, using standard capabilities in NASTRAN, for 
computing the solution of general time-dependent structural acoustics problems. Although illustrated 
for the simple geometry of spherical shell scattering, there is no restriction in the approach to 
particular geometries, so that any structure which can be modeled with NASTRAN can be handled. 
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Non— dimensionol Time 


Fig. 5. Time histories of shell velocity in the z direction. 



Non— dimensional Time 


Fig. 6. Time histories of scattered pressure. 

One of the major benefits of analyzing transient fluid-structure interaction with a general- 
purpose code like NASTRAN is the ability to integrate the acoustic analysis of a structure with other 
dynamic and stability analyses. Thus the same finite element model can often be used for modal 
analysis, frequency and transient response analysis, linear shock analysis, and underwater acoustic 
analysis. In addition, many of the pre- and postprocessors developed for use with NASTRAN 
become available for acoustics as well. 
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